{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# From 2D acoustic to 2D SH finite difference modelling\n",
    "\n",
    "In the last [lesson](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/06_2D_SH_Love_wave_modelling/1_2D_SH_FD_staggered.ipynb), we derived a 2nd order space/time finite difference solution on a staggered Cartesian grid to solve the 2D isotropic elastic SH problem. In the next step we will implement this FD scheme starting from our 2D acoustic code ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FD solution of 2D isotropic SH problem\n",
    "\n",
    "The 2nd order space/time FD scheme on a staggered grid \n",
    "\n",
    "<img src=\"../images/SG_SH-Cart.png\" width=\"50%\">\n",
    "\n",
    "using explicit time-stepping with the Leapfrog method to solve the 2D SH problem can be written as \n",
    "\n",
    "\\begin{align}\n",
    "\\scriptsize \n",
    "v_y(i,j,n+1/2) &\\scriptsize= \\frac{dt}{\\rho}\\biggl\\{\\biggl(\\frac{\\partial \\sigma_{yx}}{\\partial x}\\biggr)^c(i,j,n) + \\biggl(\\frac{\\partial \\sigma_{yz}}{\\partial z}\\biggl)^c(i,j,n) + f_y(i,j,n)\\biggr\\} + v_y(i,j,n-1/2), \\notag\\\\\n",
    "\\scriptsize\\sigma_{yx}(i+1/2,j,n+1) &\\scriptsize= dt\\; \\mu_x \\biggl(\\frac{\\partial v_{y}}{\\partial x}\\biggr)^c(i+1/2,j,n+1/2) + \\sigma_{yx}(i+1/2,j,n),\\notag\\\\\n",
    "\\scriptsize\\sigma_{yz}(i,j+1/2,n+1) &\\scriptsize= dt\\; \\mu_z \\biggl(\\frac{\\partial v_{y}}{\\partial z}\\biggr)^c(i,j+1/2,n+1/2) + \\sigma_{yz}(i,j+1/2,n), \\notag\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with the spatial derivatives on the RHS approximated by\n",
    "\n",
    "\\begin{equation}\n",
    "\\scriptsize \\biggl(\\frac{\\partial \\sigma_{yx}}{\\partial x}\\biggr)^c(i,j) \\approx \\frac{\\sigma_{yx}(i+1/2,j) - \\sigma_{yx}(i-1/2,j)}{dx},\\; \n",
    "\\biggl(\\frac{\\partial \\sigma_{yz}}{\\partial z}\\biggr)^c(i,j) \\approx \\frac{\\sigma_{yz}(i,j+1/2) - \\sigma_{yz}(i,j-1/2)}{dz},\\notag\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\scriptsize \\biggl(\\frac{\\partial v_{y}}{\\partial x}\\biggr)^c(i+1/2,j) \\approx \\frac{v_y(i+1,j) - v_{y}(i,j)}{dx},\\notag \n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\scriptsize \\biggl(\\frac{\\partial v_{y}}{\\partial z}\\biggr)^c(i,j+1/2) \\approx \\scriptsize \\frac{v_y(i,j+1) - v_{y}(i,j)}{dz}.\\notag \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the harmonically averaged shear moduli:\n",
    "\n",
    "\\begin{align}\n",
    "\\mu_{x} &= \\mu(i+1/2,j) = 2 \\biggl[\\biggl(\\mu^{-1}_{i,j}+\\mu^{-1}_{i+1,j}\\biggr)\\biggr]^{-1} \\notag\\\\\n",
    "\\mu_{z} &= \\mu(i,j+1/2) = 2 \\biggl[\\biggl(\\mu^{-1}_{i,j}+\\mu^{-1}_{i,j+1}\\biggr)\\biggr]^{-1} \\notag\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initial and boundary conditions\n",
    "\n",
    "Because we have analytical solutions for wave propagation in homogeneous media, we should test our first code implementation for a similar medium, by setting density $\\rho$ and shear modulus $\\mu$ to constant values $\\rho_0,\\; \\mu_0$\n",
    "\n",
    "\\begin{align}\n",
    "\\rho(i,j) &= \\rho_0 \\notag \\\\\n",
    "\\mu(i,j) &= \\mu_0 = \\rho_0 V_{s0}^2\\notag\n",
    "\\end{align}\n",
    "\n",
    "at each spatial grid point $i = 0, 1, 2, ..., nx$; $j = 0, 1, 2, ..., nz$, in order to compare the numerical with the analytical solution. For a complete description of the problem we also have to define initial and boundary conditions. The **initial condition** is \n",
    "\n",
    "\\begin{equation}\n",
    "v_y(i,j,0) = \\sigma_{yx}(i+1/2,j,0) = \\sigma_{yz}(i,j+1/2,0) = 0, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "so the modelling starts with zero particel velocity and shear stress amplitudes at each spatial grid point. As **boundary conditions**, we assume \n",
    "\n",
    "\\begin{align}\n",
    "v_y(0,j,n) &= \\sigma_{yx}(1/2,j,n) = \\sigma_{yz}(0,j+1/2,n) = 0, \\nonumber\\\\\n",
    "v_y(nx,j,n) &= \\sigma_{yx}(nx+1/2,j,n) = \\sigma_{yz}(nx,j+1/2,n) = 0, \\nonumber\\\\\n",
    "v_y(i,0,n) &= \\sigma_{yx}(i+1/2,0,n) = \\sigma_{yz}(i,1/2,n) = 0, \\nonumber\\\\\n",
    "v_y(i,nz,n) &= \\sigma_{yx}(i+1/2,nz,n) = \\sigma_{yz}(i,nz+1/2,n) = 0, \\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "for all time steps n. This **Dirichlet boundary condition**, leads to artifical boundary reflections which would obviously not describe a homogeneous medium. For now, we simply extend the model, so that boundary reflections are not recorded at the receiver positions.\n",
    "\n",
    "Let's implement it ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import Libraries \n",
    "# ----------------\n",
    "import numpy as np\n",
    "from numba import jit\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "\n",
    "# Ignore Warning Messages\n",
    "# -----------------------\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Definition of modelling parameters\n",
    "# ----------------------------------\n",
    "xmax = 500.0 # maximum spatial extension of the 1D model in x-direction (m)\n",
    "zmax = xmax  # maximum spatial extension of the 1D model in z-direction(m)\n",
    "dx   = 1.0   # grid point distance in x-direction\n",
    "dz   = dx    # grid point distance in z-direction\n",
    "\n",
    "tmax = 0.502   # maximum recording time of the seismogram (s)\n",
    "dt   = 0.0010  # time step\n",
    "\n",
    "vp0  = 580.   # P-wave speed in medium (m/s)\n",
    "\n",
    "# acquisition geometry\n",
    "xr = 330.0 # x-receiver position (m)\n",
    "zr = xr    # z-receiver position (m)\n",
    "\n",
    "xsrc = 250.0 # x-source position (m)\n",
    "zsrc = 250.0 # z-source position (m)\n",
    "\n",
    "f0   = 40. # dominant frequency of the source (Hz)\n",
    "t0   = 4. / f0 # source time shift (s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of 2D finite difference with analytical solution for homogeneous Vs model\n",
    "\n",
    "In a previous [exercise](https://danielkoehnsite.files.wordpress.com/2018/04/ex_02_tew21.pdf) you proved that the analytical solutions for the homogeneous 2D acoustic and 2D SH problem, beside a density factor $1/\\rho_0$, are actual identical . In the function below we solve the homogeneous 2D SH problem by centered 2nd order spatial/temporal difference operators and compare the numerical results with the analytical solution: \n",
    "\n",
    "\\begin{equation}\n",
    "u_{y,analy}(x,z,t) = G_{2D} * S \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "with the 2D Green's function:\n",
    "\n",
    "\\begin{equation}\n",
    "G_{2D}(x,z,t) = \\dfrac{1}{2\\pi \\rho_0 V_{s0}^2}\\dfrac{H\\biggl((t-t_s)-\\dfrac{|r|}{V_{s0}}\\biggr)}{\\sqrt{(t-t_s)^2-\\dfrac{r^2}{V_{s0}^2}}}, \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "where $H$ denotes the Heaviside function, $r = \\sqrt{(x-x_s)^2+(z-z_s)^2}$ the source-receiver distance (offset) and $S$ the source wavelet. Keep in mind that the stress-velocity code computes the **particle velocities** $\\mathbf{v_{y,analy}}$, while the analytical solution is expressed in terms of the **displacement** $\\mathbf{u_{y,analy}}$. Therefore, we have to take the first derivative of the analytical solution, before comparing the numerical with the analytical solution:\n",
    "\n",
    "\\begin{equation}\n",
    "v_{y,analy}(x,z,t) = \\frac{\\partial u_{y,analy}}{\\partial t} \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "To implement the 2D SH code, we first introduce functions to update the particle velocity $v_y$ ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# REPLACE SPATIAL ACOUSTIC FD BY VY UPDATE!\n",
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def update_d2px_d2pz(p, dx, dz, nx, nz, d2px, d2pz):\n",
    "    \n",
    "    for i in range(1, nx - 1):\n",
    "        for j in range(1, nz - 1):\n",
    "                \n",
    "            d2px[i,j] = (p[i + 1,j] - 2 * p[i,j] + p[i - 1,j]) / dx**2                \n",
    "            d2pz[i,j] = (p[i,j + 1] - 2 * p[i,j] + p[i,j - 1]) / dz**2\n",
    "        \n",
    "    return d2px, d2pz "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... update the shear stress components $\\sigma_{yx}$ and $\\sigma_{yz}$ ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ADD STRESS UPDATE HERE!    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and harmonically averaging the shear moduli ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ADD HARMONICALLY AVERAGING OF SHEAR MODULI HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can modify the main FD code ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     42
    ]
   },
   "outputs": [],
   "source": [
    "# 2D Wave Propagation (Finite Difference Solution) \n",
    "# ------------------------------------------------\n",
    "def FD_2D_acoustic(dt,dx,dz,f0,xsrc,zsrc):\n",
    "        \n",
    "    nx = (int)(xmax/dx) # number of grid points in x-direction\n",
    "    print('nx = ',nx)\n",
    "    \n",
    "    nz = (int)(zmax/dz) # number of grid points in x-direction\n",
    "    print('nz = ',nz)\n",
    "            \n",
    "    nt = (int)(tmax/dt) # maximum number of time steps            \n",
    "    print('nt = ',nt)\n",
    "    \n",
    "    ir = (int)(xr/dx)      # receiver location in grid in x-direction    \n",
    "    jr = (int)(zr/dz)      # receiver location in grid in z-direction\n",
    "    \n",
    "    isrc = (int)(xsrc/dx)  # source location in grid in x-direction\n",
    "    jsrc = (int)(zsrc/dz)  # source location in grid in x-direction\n",
    "\n",
    "    # Source time function (Gaussian)\n",
    "    # -------------------------------\n",
    "    src  = np.zeros(nt + 1)\n",
    "    time = np.linspace(0 * dt, nt * dt, nt)\n",
    "\n",
    "    # 1st derivative of a Gaussian\n",
    "    src  = -2. * (time - t0) * (f0 ** 2) * (np.exp(- (f0 ** 2) * (time - t0) ** 2))\n",
    "\n",
    "    # Analytical solution\n",
    "    # -------------------\n",
    "    G    = time * 0.\n",
    "\n",
    "    # Initialize coordinates\n",
    "    # ----------------------\n",
    "    x    = np.arange(nx)\n",
    "    x    = x * dx       # coordinates in x-direction (m)\n",
    "\n",
    "    z    = np.arange(nz)\n",
    "    z    = z * dz       # coordinates in z-direction (m)\n",
    "    \n",
    "    # calculate source-receiver distance\n",
    "    r = np.sqrt((x[ir] - x[isrc])**2 + (z[jr] - z[jsrc])**2)\n",
    "    \n",
    "    for it in range(nt): # Calculate Green's function (Heaviside function)\n",
    "        if (time[it] - r / vp0) >= 0:\n",
    "            G[it] = 1. / (2 * np.pi * vp0**2) * (1. / np.sqrt(time[it]**2 - (r/vp0)**2))\n",
    "    Gc   = np.convolve(G, src * dt)\n",
    "    Gc   = Gc[0:nt]\n",
    "    lim  = Gc.max() # get limit value from the maximum amplitude\n",
    "    \n",
    "    # Initialize empty pressure arrays\n",
    "    # --------------------------------\n",
    "    p    = np.zeros((nx,nz)) # p at time n (now)\n",
    "    pold = np.zeros((nx,nz)) # p at time n-1 (past)\n",
    "    pnew = np.zeros((nx,nz)) # p at time n+1 (present)\n",
    "    d2px = np.zeros((nx,nz)) # 2nd spatial x-derivative of p\n",
    "    d2pz = np.zeros((nx,nz)) # 2nd spatial z-derivative of p\n",
    "\n",
    "    # Initialize model (assume homogeneous model)\n",
    "    # -------------------------------------------\n",
    "    vp    = np.zeros((nx,nz))\n",
    "    vp    = vp + vp0       # initialize wave velocity in model\n",
    "\n",
    "    # Initialize empty seismogram\n",
    "    # ---------------------------\n",
    "    seis = np.zeros(nt) \n",
    "    \n",
    "    # Calculate Partial Derivatives\n",
    "    # -----------------------------\n",
    "    for it in range(nt):\n",
    "    \n",
    "        # FD approximation of spatial derivative by 3 point operator\n",
    "        d2px, d2pz = update_d2px_d2pz(p, dx, dz, nx, nz, d2px, d2pz)\n",
    "\n",
    "        # Time Extrapolation\n",
    "        # ------------------\n",
    "        pnew = 2 * p - pold + vp ** 2 * dt ** 2 * (d2px + d2pz)\n",
    "\n",
    "        # Add Source Term at isrc\n",
    "        # -----------------------\n",
    "        # Absolute pressure w.r.t analytical solution\n",
    "        pnew[isrc,jsrc] = pnew[isrc,jsrc] + src[it] / (dx * dz) * dt ** 2\n",
    "                \n",
    "        # Remap Time Levels\n",
    "        # -----------------\n",
    "        pold, p = p, pnew\n",
    "    \n",
    "        # Output of Seismogram\n",
    "        # -----------------\n",
    "        seis[it] = p[ir,jr]   \n",
    "        \n",
    "    # Compare FD Seismogram with analytical solution\n",
    "    # ---------------------------------------------- \n",
    "    # Define figure size\n",
    "    rcParams['figure.figsize'] = 12, 5\n",
    "    plt.plot(time, seis, 'b-',lw=3,label=\"FD solution\") # plot FD seismogram\n",
    "    Analy_seis = plt.plot(time,Gc,'r--',lw=3,label=\"Analytical solution\") # plot analytical solution\n",
    "    plt.xlim(time[0], time[-1])\n",
    "    #plt.ylim(-lim, lim)\n",
    "    plt.title('Seismogram')\n",
    "    plt.xlabel('Time (s)')\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    plt.show()            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "dx   = 1.0   # grid point distance in x-direction (m)\n",
    "dz   = dx    # grid point distance in z-direction (m)\n",
    "\n",
    "dt = 0.001   # time step (s)\n",
    "\n",
    "FD_2D_acoustic(dt,dx,dz,f0,xsrc,zsrc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "* Implementation of the 2D isotropic elastic SH problem based on our 2D acoustic FD code \n",
    "* Comparison with analytical solution for homogeneous Vs model"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
